Importing additional spatial data

library(sf)
library(mapview)
library(dplyr)

Mined Units

The mined units can be read in directly with st_read.

mined <- st_read('https://ca.dep.state.fl.us/arcgis/rest/services/OpenData/MMP_MINEDUNITS/MapServer/0/query?outFields=*&where=1%3D1&f=geojson')
Reading layer `OGRGeoJSON' from data source 
  `https://ca.dep.state.fl.us/arcgis/rest/services/OpenData/MMP_MINEDUNITS/MapServer/0/query?outFields=*&where=1%3D1&f=geojson' 
  using driver `GeoJSON'
Simple feature collection with 371 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -82.95369 ymin: 27.46063 xmax: -81.69582 ymax: 30.51925
Geodetic CRS:  WGS 84
mapview(mined)

Brownfield sites

Download the zipped kml file to a temporary directory.

# url with zipped kml
urlin <- 'https://ordsext.epa.gov/FLA/www3/acres_frs.kmz'

# download file
tmp1 <- tempfile(fileext = ".kmz")
download.file(url = urlin, destfile = tmp1, method = 'curl')

Unzip the kmz file.

tmp2 <- tempdir()
unzip(tmp1, exdir = tmp2)

Get the name of the kml file to read.

lyr <- unzip(tmp1, list = T)$Name
fl <- paste(c(tmp2, lyr), collapse = "\\")
fl <- gsub('\\\\', '/', fl)

Read the kml file with st_read and drop Z dimension with st_zm. Here, the Tampa sites are loaded. You can view all possible locations in the kml file with st_layers.

dat <- st_read(fl, layer = 'TAMPA') %>% 
  st_zm()
Reading layer `TAMPA' from data source `/tmp/RtmpqdPSoL/ACRES_FRS.KML' using driver `LIBKML'
Simple feature collection with 103 features and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: -82.51857 ymin: 27.87108 xmax: -82.3533 ymax: 28.07826
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
mapview(dat)

To import all layers in the kml file, identify the layer names and loop through them to add to a single object. This takes several hours. The data are saved as an .Rdata object and .csv file for later use. The data include only the site name and location in decimal degrees.

# layer names
alllyr <- st_layers(fl)$name

strt <- Sys.time()
out <- NULL
for(i in alllyr){
  
  # counter
  cat(i, which(i == alllyr), 'of', length(alllyr), '\n')
  print(Sys.time() - strt)
  
  # import each layer
  dat <- st_read(fl, i, quiet = T)[, c('Name')]
  
  # append to same object
  out <- rbind(out, dat)
  
}

# save as RData object
allbfld <- out %>% st_zm()
save(allbfld, file = 'data/allbfld.RData', compress = 'xz')

# save as csv
allbfldcsv <- allbfld %>% 
  mutate(
    lon = st_coordinates(.)[,1], 
    lat = st_coordinates(.)[,2]
  ) %>% 
  st_set_geometry(NULL)
write.csv(allbfldcsv, 'data/allbfldcsv.csv', row.names = F)

View all brownfield sites.

load(file = 'data/allbfld.RData')
mapview(allbfld, legend = F, col.regions = 'grey')

Unlink the temporary files to delete them when you are finished.

unlink(tmp1, recursive = TRUE)
unlink(fl, recursive = TRUE)

Census tracts

This workflow demonstrates how to download census tracts for the United States.

Download the zip file to a temporary directory and unzip to a second temporary directory.

urlin <- 'https://static-data-screeningtool.geoplatform.gov/data-versions/1.0/data/score/downloadable/1.0-shapefile-codebook.zip'

tmp1 <- tempfile(fileext = ".zip")
download.file(url = urlin, destfile = tmp1)

tmp2 <- tempdir()
unzip(tmp1, exdir = tmp2)

The usa.zip file was included in the previously zipped file. Unzip usa.zip file.

zip1 <- list.files(tmp2, 'usa\\.zip', full.names = T)
unzip(zip1, exdir = tmp2)

Get the file path for usa.shp and import with the sf package.

fl <- list.files(tmp2, '\\.shp', full.names = T)
tracts <- st_read(fl)
Reading layer `usa' from data source `/tmp/RtmpqdPSoL/usa.shp' using driver `ESRI Shapefile'
Simple feature collection with 74134 features and 123 fields (with 367 geometries empty)
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -179.2311 ymin: -14.60181 xmax: 179.8597 ymax: 71.44106
Geodetic CRS:  WGS 84

Clip to the Tampa Bay watershed boundaries using the tbshed spatial object. The CRS is the same, so there is no need to transform.

load(file = 'data/tbshed.RData')

tbtracts <- tracts %>% 
  st_intersection(tbshed)

save(tbtracts, file = 'data/tbtracts.RData')

View the data.

mapview(tbtracts)

Remove temporary files.

unlink(tmp1, recursive = TRUE)
unlink(zip1)
fls <- list.files(tmp2, gsub('\\.shp$', '', basename(fl)), full.names = T)
file.remove(fls)
[1] TRUE TRUE TRUE TRUE TRUE